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Abstract 

The Beauchemin model is a simple particle-based description of stochastic lymphocyte migration in tissue, which 
has been successfully applied to studying immunological questions. In addition to being easy to implement, the 
model is also to a large extent mathematically tractable. This article provides a comprehensive overview of both 
existing and new analytical results on the Beauchemin model within a common mathematical framework. 
Specifically, we derive the motility coefficient, the mean square displacement, and the confinement ratio, and 
discuss four different methods for simulating biased migration of pre-defined speed. The results provide new 
insight into published studies and a reference point for future research based on this simple and popular 
lymphocyte migration model. 



Background 

A unique property of the immune system is that it mainly 
consists of constantly moving cells. Technological progress 
has facilitated the study of lymphocyte migration at 
increasingly higher resolutions, culminating in the applica- 
tion of two-photon microscopy to tracking single lympho- 
cytes in the living animal [1-3]. Accompanying this 
progress, mathematical and computational models of lym- 
phocyte migration have been developed and applied to 
both qualitative and quantitative questions. These models 
can be roughly categorized as follows according to their 
level of detail: Whole-population models, usually formu- 
lated using ordinary or partial differential equations, do 
not distinguish between individual cells, but treat cell sub- 
populations as continuous quantities. Particle-based mod- 
els do consider each cell individually, but treat cells as 
freely moving particles without mass, volume, and shape. 
Finally, whole-cell models such as the Cellular Potts 
Model [4] explicitly represent the cell and its interaction 
with the environment. In the present paper, we analyze a 
particle-based model that was introduced by Beauchemin, 
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Dixit, and Perelson [5] . Throughout the paper, we refer to 
this model as the "Beauchemin Model". 

Particle-based approaches such as the Beauchemin 
model are used for studying questions where a whole- 
population approach does not provide sufficient informa- 
tion, but a whole-cell model is not necessary, not feasible 
for computational reasons or would require too many 
assumptions on unknown parameters. For instance, 
Grigorova et al. [6] and ourselves [7] recently used the 
Beauchemin model to study the transit of T cells through 
a lymph node in the absence of antigen, and we used it to 
determine the amount of directional bias that could be 
detected in random T cell migration using contemporary 
two-photon imaging techniques [7] . 

The Beauchemin model is defined by three parameters - 
a speed Vf ree , a time interval tf ree , and another time interval 
£ P ause- The model describes a particle moving freely in a 
three-dimensional space according to the following rules: 
The particle starts at a fixed position and turns in a ran- 
dom direction (more precisely, a direction is sampled from 
a uniform distribution on the unit sphere). The turn takes 
a fixed time t pause , during which the cell does not move. 
This reflects the time it takes to "displace the internal 
structure which brings about motion" [4] . Then, the cell 
moves in the chosen direction with velocity v free for a fixed 
time if ree , after which it stops and restarts the process. 
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The main distinguishing characteristic of the model is 
the introduction of the pause phase. In this sense, it can 
be seen as a generalization of the "ideal chain" model of 
polymer physics [8], in which a chain of rigid rods of 
fixed length t is considered that are freely jointed to each 
other. This is mathematically equivalent to the Beauche- 
min model if we set v free • tfree = I and f pause = 0. In par- 
ticular, analytical results e.g. on the mean square 
displacement of the ideal chain [8] are directly transfer- 
able to this case. As we shall see later on in this paper, 
for some quantitative properties even the more general 
case with £ pause 5 0 can be reduced to the ideal chain 
model. 

This paper is structured as follows: In the upcoming sec- 
tion, we begin by formulating a mathematical framework 
for our analysis, which explicates sufficiently but not 
overly general notions of random walks and Brownian 
motion. We deliberately include a fairly large amount of 
detail, including proofs, in order to present our results in a 
self-contained fashion. Thereafter, we present analytical 
results about random motion in the Beauchmin model, 
and discuss four different ways for simulating biased 
migration with pre-defined speed. We wrap up with con- 
cluding remarks. 

We point out that around half of the results discussed 
in this paper originate from our recent study that used 
the Beauchemin model [7]. We chose to integrate these 
results into the present paper in rewritten and extended 
form because they were formerly only published in con- 
densed form as supporting online material. In this man- 
ner, we hope to provide a coherent account of all 
available analytical results, which we anticipate to serve 
as a useful reference for future users of the model. We 
mark results that have appeared earlier in our previous 
study or elsewhere (e.g. Proposion 10) by giving the 
appropriate reference in proposition or theorem state- 
ments, while propositions or theorems stated without 
such a reference (e.g. Proposion 12) have, to our knowl- 
edge, not been published before. 

Mathematical foundations 

In this section, we outline the basic mathematical frame- 
work that we use to describe stochastic cell motion and, 
later on, the Beauchemin model. A similar, but slightly 
less general framework (based on random walk on a lat- 
tice) is used in Berg's textbook Random Walks in Biology 
[9], while a much more general account of random walks 
than given here is usually presented in probability text- 
books (e.g., [10]). At the end of this section, we will outline 
the connection of this random walk to the corresponding 
macroscopic partial differential equation model, which is 
the convection-diffusion-equation. 

First of all, we introduce some notational conventions. 
"We denote the expectation of a random variable f by E 



its variance by Var(f) = E [(£ - E[f]) 2 ] and its prob- 
ability density function by f%. Elementary random vari- 
ables are written in lowercase Greek letters (usually tf ; ) 
while derived random variables that are functions of ele- 
mentary random variables are written in capital Latin let- 
ters. All random variables are assumed to be continuous 
and real-valued. Mean and variance of random variables 
will also be written as \i and a 2 , if the referenced variable 
is clear from context. 

Observed tracks and associated measures 

The Beauchemin model is usually validated against, or 
used to generate, cell tracking data as it would be 
recorded in a two-photon microscopy experiment. As 
we will see later, there can be subtle differences between 
the real underlying state and the observed state of a cell. 
We thus make this difference explicit in our analysis by 
means of the following definitions for cell tracking 
terminology. 

Definition 1 (Cell track). A d-dimensional cell track is 

a finite sequence 

(x(l),t(l)),..., (x(fe),t(fe)) 

consisting of positions x(l) , . . . , x{k) e M. d and increas- 
ing time points t(l), . . . , t(fe) e R. 

The time indices are written in brackets to indicate 
that the positions are discrete-time observations of an 
underlying continuous motion process. For population- 
based cell track analysis, it is common to align a set of 
tracks to a common starting point. 

Definition 2 (Aligned cell track sets). For a cell track 
T = {x(l), t(l)) , (x(k), t{k)), the track that results 
from subtracting the first element from all elements in 
the sequence, i.e. 

T 0 = (0,0),(x(2)-x(l),t(2)-t(l)) , (x(k)-x(\),t(k) - t(l)) 

is called the zero-aligned version of T . 

Throughout the paper, we will assume that all cell 
track sets are zero-aligned. 

Definition 3 (Mean displacement and mean square 
displacement). Let X\(f) , x n {t) be the positions of n 
zero-aligned cell tracks at some fixed time t. Then the 
mean displacement D{f) for S is defined by 

DW = + ■■■ +]M0N 

n 

where \ \ ■ \ \ is the d-dimensional euclidian norm or 
vector length, and the mean square displacement D 2 (t) is 
defined by 

ra,., Il*i(0ll 2 + ••• +IW0II 2 
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Beauchemin et al. [5] estimated the three parameters 
Vfree> ^free and £ paU se of their model by fitting it to mean 
displacement data: For each combination of parameters, 
10 6 simulations of the model were run, and the mean dis- 
placement data was recorded and averaged over all simu- 
lations. The resulting mean displacement curve was 
compared to real data from several publications. The 
parameter triplets were ranked according to their 
squared sum of errors to the real data (see Table 1). It 
was found that the values f free = 2.0 min, f paU se = 0-5 min 
and v free = 18.8 ^m/min gave the best agreement to the 
real data. These values are biologically reasonable and 
similar to those estimated by Miller et al. [1]. 

Microscopic random walk with bias 

In this section, we will define the notion of a one-dimen- 
sional random walk used in this paper and derive its rela- 
tion to a one-dimensional convection-diffusion process. 
These concepts can be easily generalized to higher dimen- 
sions under the assumption that the dimensions are pair- 
wise uncorrelated, while stochastic independence in the 
strict sense is not necessary. This assumption holds for the 
particular random walks we analyze in this paper. 

Definition 4 (One-dimensional random walk [9,11]). 
Let i gNi be independent, identically distributed real- 
valued random variables with mean n and variance a 2 . 
Then the sequence (S t ), i e N, defined by 

i=i 

is called a random walk with step bias \i and step var- 
iance a 2 . If n =0, (S t ) is called an unbiased random 
walk. 

Intuitively, the sequence (S t ) is the random walk tra- 
jectory and the are the individual steps. The above 
definition is a slightly generalized version of the com- 
monly used entirely discrete random walk with fixed 
steps to the left and the right: We assume that each 



Table 1 Different parameters can lead to similar motility 
coefficients. 



Rank 


^pause 


fyree 


v free 


M 


SSR 




min 


min 


fjm/min 


f/m 2 /min 




1st 


0.5 


2.0 


18.8 


94.25 


1.00 


2nd 


0.5 


2.5 


16.6 


95.68 


1.01 


3rd 


0.25 


2.0 


17.6 


91.78 


1.02 


4th 


1.5 


1.5 


26.1 


92.89 


1.03 


5th 


0.75 


1.5 


23.8 


94.41 


1.03 



Parameter triplets of the 3D cell random walk model that were found to best 
fit in vivo data by Beauchemin et al. [5] along with their corresponding 
motility coefficients as per equation (2). Beauchemin et al. calculated sums of 
squared residuals (SSR) to a set of in vivo mean displacement data from four 
publications [24,25,1,26]. The SSR values shown here are normalized to the 
best ranking triplet. 



step is sampled from a (possibly continuous) random 
variable with finite first and second moments. The fol- 
lowing corollary provides an important characteristic of 
this random walk: 

Corollary 5 (Mean and variance of a random walk). 
For a random walk (S t ) as given by Definition 4, 

E[S t ] = f/x and Var(S t ) = to 2 . 

Proof. Since the cf, are independent and thus uncorre- 
lated, this is a direct consequence of the linearity of var- 
iance and expectation. □ 

If the random walk is unbiased, Var(S t ) can be geome- 
trically interpreted as the expected mean square squared 
displacement of the particles at time t (Definition 3). 
This leads to an important characteristic of the unbiased 
random walk: In a particle ensemble, the particles will 
go nowhere on average, and their expected mean square 
displacement is linear in time. This characterization is 
often used in the immunological literature [2] to convey 
the difference between a random walk versus directed 
migration, where the root mean square displacement is 
linear in time (Figure 1). 

For sufficiently long observations, we do not only 
know the variance and mean of a random walk, but can 
furthermore obtain a good approximation of the entire 
particle distribution from the central limit theorem. 

Theorem 6 (Central limit theorem [9,10]). Let (S t ) be 
a random walk as given in Definition 4. Then as t —> °°, 
the distribution of the random variable 

7 S t -tn 
A - -p — 

converges to a normal distribution with mean 0 and 
variance a 2 . 

Thus, after a sufficiently long time, the distribution of 
the random walk (S t ) at fixed times is approximately 
normal. Of course, the terms "sufficiently long" and 
"approximately normal" are rather vague. In the context 
of intravital two-photon imaging, which is limited to a 
finite, often axially thin imaging region, we can usually 
not safely assume that we are in the range of validity of 
the normal distribution approximation: We are obser- 
ving only a few "steps" (prolonged periods of rather per- 
sistent motion followed by short pauses) of a random 
walk. Beltman et al. [3] and ourselves [7] have discussed 
the implications of this technical limitation for detecting 
directional bias and estimating the motility coefficient. 

We now generalize the random walk to higher 
dimensions. 

Definition 7 (Multi-dimensional random walk [7]). Let 

d €N- Suppose (f/ 1 ', . . . , §/ d ')> i e N are independent, 
identically distributed ^-valued random vectors 
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directed 
random 





time 

Figure 1 Directed vs. random motion. Traveled distance (directed movement) and root mean square displacement (random movement) of 
directed and random motion plotted as functions of time (left) and square root of time (right). On the left hand side, directed motion yields a 
straight line, with a slope given by its speed. On the right hand side, random motion yields a straight line, with a slope given by its diffusion 
coefficient (also called motitity coefficient in this paper). 



with ^\ ... , ^ pairwise uncorrelated, 

t'ft' 1 '] = ^ £[#'] = ^ and Var{^ ] ) = ...= Var{^ d) ) = a 2 - 

Then the sequence {S t ), t e N, defined by 



(i) 



)-Efc ll) r ] ) 



i=i 



is called a d-dimensional random walk with step bias 
{fA m , // d) ) and step variance a 2 . If(i m = ... = ^ <d) = 0, 
then (Sf) is called an unbiased d-dimensional random 
walk. 

The following theorem gives the limit distribution of 
suitably scaled multi-dimensional random walks within 
our framework. 

Theorem 8 (Multi-dimensional central limit theorem 
[7]). Let (S t ) be a d-dimensional random walk as given 
in Definition 7. Then as t — » °°, the distribution of the 
random vector 



(S. 



(1) 



£M (1) S 



id) 



41 



converges to a d-dimensional normal distribution with 
zero means and covariance matrix 

E = a 2 E 

Proof. Let a\, . . . , ad € R. According to the Cramer- 
Wold Theorem, it suffices to show that 



[a^SY* - t/z«) + ... + a d {Sf ] - t/iW))/yi converges 
in distribution to a normal random variable with zero 



+ a 2 A )a 2 (see [12]). Note that 



are 



mean and variance (a 2 

the random variables a x ^ (1) + ... + a&\~', te 
independent, identically distributed and have mean 
fliji/ 1 ' + ... + adft (d) . Furthermore, since f^ 1 ', . . . , ^ are 
pairwise uncorrelated, 

Var^f/ 1 ' a4 d) ) = {a 2 1+ ... +a 2 )a 2 . 



Hence, the result follows from the central limit theo- 
rem for one-dimensional random walks. D 

The central limit argument shows that even complex 
cell migration processes might have a simple description 
when observed for a sufficiently long time: The diffusion 
coefficient is a function of all "microscopic" cell motility 
parameters such as persistence time, speed, turning 
angle distributions and so on. In the scaling limit, this 
single coefficient describes an unbiased random walk 
completely. Therefore, it is more common to describe 
random walk processes on such timescales in terms of a 
differential equation instead of a stochastic process. This 
can be done using a scale transition, as we will explain 
below. 

The convection-diffusion equation 

The random walk (S t ) introduced in Definition 4 is a 
plausible microscopic model for the position of particles. 
However, if the positions are observed on large time 
scales (e.g., only every 20th state of the random walk 
can be observed due to a limited temporal resolution), 
then a macroscopic continuous model can be used to 
describe the migration. This continuous model no 
longer depends on the detailed characteristics of the 
microscopic model such as the turning angle distribu- 
tion or persistence length (i.e., the precise shape of 
the £). 

Mathematically, we construct the continuous model as 
follows: For n e N let {Z n (t)) with t > 0 be given by 



1 



|mj 



Z n [t) = -j= J] ,- n) + til. 



(1) 



We write the time index t in brackets in order to 
emphasize that the processes are given in continuous 
time. The process (Z„ (t)) is obtained by observing n steps 
of the random walk (S t ) in every time interval (t - 1, t] 
with t € N. It is scaled such that E [Z n (t)] = E [S t ] = tfi 
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and Var(Z„(f)) = Var{S t ) = to 2 for every t e N. Between 
two successive steps, the process is constant. 

Let (Z{t)), t > 0, denote the limiting process obtained 
as n, the number of steps per time interval, tends to 
infinity. According to the Functional central limit theo- 
rem (sometimes also referred to as Donsker's Theorem, 
see [12]), (Z(t)) has the following two properties: 

♦ Any finite-dimensional distribution of (Z(t)) is 
Gaussian. 

♦ E[Z(£)] = tfi and Cov(Z(s), Z{t)) = mm(s, fjcr 2 for 5, 
t > 0. 



common in the immunological literature since Miller et 
al.'s early work [1]. 

Similarly as for the discrete model, the continuous 
model can be generalized to multiple dimensions. In 
this case, (Z„ (t)) is obtained by observing n steps of a 
^-dimensional random walk (S t ) in every time interval 
of unit length. For n tending to infinity, the Functional 
central limit theorem yields convergence to a d-dimen- 
sional Brownian motion with independent components 
and convection coefficients C (1) = pi w , C (rf) = ^ d K 
The motility coefficient is equal to M = do 2 12, where d 
is the number of observed dimensions. 



A process with these properties is called Brownian 
motion with drift ft. If t>0, then Z{t) has a probability 
density <p{x, t) given by 



</>(x, i) = 



1 



V 'Into 



exp 



(x — t[iy 
2ta 2 



for x e R. 

As one can easily verify, this probability density func- 
tion solves the following well-known differential equa- 
tion for the initial condition (p(x, 0) = S(x), where d 
denotes the Dirac delta distribution. 

Theorem 9 (Convection-diffusion equation [11,13]). 
Let C = p and M = <r 2 /2. Then 



at 



<P{x, t) = 



d 3 2 
-C— +M — 

dx dx z 



(p(x, i). 



Taken the initial condition given above, the solution 
of the convection-diffusion-equation exists and is unique 
[13]. 

As the sum of independent Brownian motions is again 
a Brownian motion, the convection-diffusion equation 
describes the dynamics of single diffusing particles as 
well as of ensembles of particles carrying out motions 
independently: Assuming the ensemble is large enough 
such that stochastic perturbations can be ignored, the 
convection-diffusion-equation describes the evolution of 
the particle concentration over time. This scaling argu- 
ment justifies the use of the convection-diffusion equa- 
tion for simulating large cell populations on timescales 
that are substantially longer than the rhythmicity of the 
cell migration process. For example, we used the con- 
vection-diffusion equation to simulate the transit of 
T cells through the lymph node, which takes approxi- 
mately half a day [7]. 

Note that in the context of molecular motion, the 
quantity M is usually called diffusion coefficient. 
Throughout this paper we will use the term motility 
coefficient to indicate that we are talking about particles 
that are substantially larger than molecules; this is 



Summary 

As a microscopic model for stochastic motility of parti- 
cles, we have first introduced one-dimensional random 
walks. According to this model, particles move by taking 
independent and identically distributed steps at equidi- 
stant discrete time points. Therefore, just like in the 
Beauchemin model, our framework does not consider 
persistence to last longer than the duration of a single 
step; the model is fully characterized by the absolute 
time between successive steps and the distribution of 
steps. Important characteristics of the step distribution 
are the step bias n and the step variance a 2 . The ran- 
dom walk is called unbiased if \i = 0, and biased, other- 
wise. If we observe many particles carrying out 
simultaneous unbiased random walks, the particles will 
go nowhere on the average, but they will spread out 
increasingly, their distribution being approximately nor- 
mal. The model can be easily generalized to d dimen- 
sions if we assume that the steps are dimension-wise 
uncorrelated (not necessarily statistically independent). 
Essentially, we can then treat the d dimensions like d 
separate one-dimensional random walks. 

On the macroscopic scale, (biased) random walks can 
be approximated by Brownian motion (with drift). In 
this case the motion of particles is described sufficiently 
well by only two parameters, namely the motility coeffi- 
cient M = a 2 /2 and the convection coefficient C = fi. 
The motility coefficient M quantifies the random com- 
ponent of a biased random walk, while the convection 
coefficient C (which can be interpreted as a velocity) 
characterizes the directed component, if there is any 
(Figure 2). Convergence to Brownian motion can also be 
obtained for much more general random walk models 
using appropriate central limit theorems. However, the 
particular random walk chosen here is general enough 
for our purpose of analyzing the Beauchemin model. 

The Beauchemin model 

In this section, we analytically derive three quantitative 
characteristics of the particular random walk defined by 
the Beauchemin model, namely the motility coefficient, 
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the mean square displacement, and the confinement 
ratio. 

Motility coefficient 

The fact that the model is three-dimensional might 
seem to complicate our analysis. However, fortunately, it 
is possible to decouple the dimensions of the model, as 
will be shown in the following. The upcoming result 
was stated by ourselves [7] and, more recently, by Dono- 
van and Lythe [14]. 

Proposition 10 (Motility Coefficient of the Beauche- 
min model [7,14]). Let i>,-, f,) e R 3 , [ e N, be indepen- 
dent random vectors distributed uniformly on a sphere 
with radius r. Set (X t , Y tl Z t ) = Yl\=i v i' Then X„ 
Y t and Z t are unbiased random walks as defined in 4, 
and as t — > °°, the distribution of (X t ,Y t ,Z t )/*/t con- 
verges to a 3-variate normal distribution with zero mean 
and covariance matrix 




r 2 /3. 



Proof. The key ingredients to the proof are 

(i) Zi, Vi, Q are uniformly distributed on [-r, r], 
(it) Cov(4 vd = Cov(4 Q = Cov(vi, Q = 0. 

In order to establish (i), note that P (f, < h), i.e., the 
probability that lies below the height h in the sphere, 
is proportional to the part of the sphere's surface area 
lying below the plane = h. The surface of a spherical 
cap is proportional to its height, and thus <f, must be 
uniformly distributed on [-r, r] (see [15]). By isotropy, 
the same arguments hold for v t and Ci- The statement 
(ii) on pairwise uncorrelatedness is obtained by (ii) is 
obtained by observing that (£,, y,) is uniformly distribu- 
ted on a circle with radius r - \z\ given that Q = Z, and 
hence Cov(f ; , y,) = 0. 



Now, (i) implies E [£,•] = E [yj = E [Q] = 0, and hence 
(X t ), (Y t ) and (Z t ) are unbiased random walks. The state- 
ment on the asymptotic distribution of {X t , Y t , Z t ) follows 
by (ii) together with Theorem 8, and the fact that a ran- 
dom variable with uniform distribution on [-r, r] has 
variance rV3. a 

Hence, particles in the Beauchemin model behave 
exactly like the superposition of three uncorrected ran- 
dom walks with uniformly distributed steps. 

Mean square displacement 

Putting together our result from the last section with 
the relationship between the motility coefficient M and 
the dimension-wise step variance a 2 of a random walk, 
we obtain the following relation between the diffusion 
coefficient and the microscopic parameters: 



M = 



(^free X tfree) 
6(tfree + tpause) 



(2) 



The corresponding diffusion coefficients for the 5 
best-fitting parameter triplets determined by Beauche- 
min et al. are given in Table 1. It turns out that all the 
best-fitting triplets have very similar motility coefficients. 
On the other hand, the microscopic motility parameters 
vary quite a lot: For example, the pause time of the 4th 
best fitting triplet is 5 times longer than that of the 3rd 
best fitting triplet. The triplets were ranked by Beauche- 
min et al. according to the sum of squared residuals 
(SSR) to the experimental data. However, the SSR of the 
5th best fitting triplet is only less than 4% higher than 
that of the best fitting triplet. Our upcoming analysis 
will provide some insight why this occurs. 

The next two propositions will yield a precise charac- 
terization of the variance - and with it the mean square 
displacement - of the position of single particles and 
ensembles of observed particles in the Beauchemin 
model. It is mathematically more convenient to analyze 
the mean square displacement rather than the mean dis- 
placement because, as we have shown above, the 




x 0 —Or x 0 

Figure 2 Position distributions for biased random walks. Expected distribution of zero-aligned cell displacements measured in one 
dimension in the case of an unbiased random walk with motility coefficient M (left) and a biased random walk with convection coefficient C 
(right) as per theorem 9. In both cases, the distribution is Gaussian with standard deviation ^/2tM- ln tne biased case, the mean of the 
distribution moves away from the origin with time. 
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dimensions can be decoupled. For mean displacement, 
this does not hold, and we would have to take into 
account the stochastical dependence between 
dimensions. 

We start out with the variance of the position of one 
particle that starts its random walk at time t = 0. For sim- 
plicity, we will derive the following equations assuming 
that one dimension is observed. Due to the independence 
of the dimensions discussed above, multi-dimensional 
versions of these results are obtained simply by multiply- 
ing the derived quantity (i.e., mean square displacement 
or confinement ratio) with the number of observed 
dimensions. 

Proposition 11 (Expected square displacement of sin- 
gle particles in the Beauchemin model). A particle starts 
at time t = 0 in position x{0) = 0 to carry out a random 
walk as defined by the Beauchemin model with fixed 
parameters t pause ,tfr ee ,Vj i . ee . Let x(t + r) denote the position 
of the particle observed in one dimension at time t + r 
where t is an integer multiple of tf ree + t pause and 0 < r < 
tfree + tpause- Then the expected square displacement of 
the particle is given by: 



E [D 2 [t+T)\ : 



Var(x{t + r)) = 2Mt + -^(max(t ■ 



tpausei 0)) 



with M as defined in equation (2). 

Proof. The term 2Mt follows directly from Proposition 
10 and the linearity of variance (see also Corollary 5). 
Thus we continue with t = 0. Then we have from Pro- 
position 10 that the position at time t free + t pauS e is 



uniformly distributed on the interval [-v bee t f[ee , v free t free ], 
having variance 2M(t free + f pause )- 

For r < £ P ause» the particle's position does not change, 
hence its variance is equal to 0. For r > £ pause , the posi- 
tion is uniformly distributed on an interval with a size 
proportional to r - t pause - Hence the variance has quad- 
ratic form y 
we obtain 



a(t - t pRUSe ) • Inserting x = t bee + f, 



pause' 



2M( tfree + tpause) = « tfree 

Inserting the definition of M from above, this results 



in 



a = 



(ffree) 



(3) 



Combining the cases r < t pause and r > t pause , we arrive 
at the proposition. □ 

A plot of the expected square displacement shows a 
series of iterated pulses that approaches the linear pro- 
gression described by the diffusion equation from below 
(see Figure 3). However, we will most likely not see 
such phasic behaviour if we take the average of several 
observed particles, since the particles will in general not 
be synchronized and the pulses that are caused by the 
free runs between turns will average out. Consequently, 
the next step in our analysis is to derive the expected 
mean square displacement for a particle ensemble. We 
do this by assigning to each observed particle a phase, 
which reflects the particle's state at the time when we 



o 
o 

I 



exact variance 
diffusion equation 



0 tc 



time 



Figure 3 The Beauchemin model and the diffusion equation. The expected square displacement (or position variance) of a single particle 
along with the approximation by the diffusion equation. 
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start our observation. Then we average over all possible 
phases. 

Proposition 12 (Empirical mean square displacement 
of the Beauchemin model). Set t run = tf ree + t pause A par- 
ticle starts at a uniformly distributed random time -t 0 , 
-trun <-to < 0 in some position to carry out a random 
walk as defined by the Beauchemin model with fixed 
parameters t pause> tf ree> Vf ree . Let x(t) e R denote the posi- 
tion of this particle observed in one dimension at time t, 
where we choose the coordinate system such that x(0) = 
0. Then the expected square displacement of the particle 
at time t is given by 



E[D 2 (t)] 



trun J-tnu 



Var(x(t - t 0 ) - x(t 0 )) dt 0 



= 2Mt - 2Mtfr ee X 



3 f — tf ree 



Proof. See Figure 4 for an illustration of how the para- 
meter 8 affects the variance of the observed particle 
position during the initial observation time. For 8 < 0, 
the cell is resting at the beginning of observation time 



x(t) = x{t + r 0 ) — x(t 0 ) 

where x(t) denotes our observed position. We are 
interested in the variance Vflr[x(t)], for which the follow- 
ing holds: 



Var[x(t)] 

= Var[x(t + r 0 ) — x(r 0 )] 



(4) 
(5) 



= Var[x(t + t 0 )] + Var[x(r 0 )] — 2Cov[x(t + to, x(r 0 )] (6) 

What we are deriving is the expectation of expression 6, 
which is a function of the random variable r 0 . Due to the 
markov property of the random walk, only the last 
observed step matters and we can take t 0 as uniformly dis- 
tributed OVer [0, ff ree + f paU se)' Thus, 



E[Var[*(t)]j . 



tfree 



/ Var[x(t+ t 0 )] - Var[x(r 0 )]dt 0 

+ tpause Jz=0 



We start by integrating the two variance terms. Some 
algebra yields 

tfree+tpause 



the cell is resting at the beginning ot observation time 1 /-ftee+Jpause 

and remains paused for time £ pause - 8. For 8 >0, the 7~~t / Var[x(t + r 0 )] + Var[r 0 ]dr 0 

, pause > tfee + tpause Jt 0 =0 



cell is moving at the beginning of observation time and 
keeps moving for time 8. Let us write Var(Ss{t)) for the 
position variance of a particle with phase shift 8. 

Denote by r 0 the time that the particle was already 
moving when we started observing it. Hence, 



( Vfree tfree) 



3 (tfree + tpause) 

= 2 M Ut + t ^j + t ^\ 



tfree \ tfree 
t + + 

3 J 3 



r5 < 0 



tpause 

<5 = 0 




0 



^free H~ 



8 > 0 




0 



^frcc tv; 



time time 
Figure 4 Effect of the phase parameter on the variance. The plots illustrate how the phase parameter d affects the variance of a single 
observed cell in the Beauchemin model. The expected variance for a (large) ensemble of observed particles is obtained by averaging over all 
values of 5. 
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where M is the motility coefficient as defined in equa- 
tion 2. For the covariance term it can be shown by 
some case distinctions that 



1 /"'free+fpause 

/ Cov[x(t+ r 0 ), r 0 \dr 0 

+ fpause Jz a =0 



tfrei 



mm(t, tftee) min([, t free )3 2( fre , 
= M min(t, t free ) + M I + — 2 + — 

tfree M flee i 



=0 for (>tfr«. 

Combining the variance and covariance terms as per 
expression 6, we obtain the claimed identity. n 

As illustrated in Figure 5, the empirical mean square 
displacement from time tf lee onwards is thus perfectly 
described by the linear equation 



E [D 2 (t)] =at + fi 
with slope 

(Vfree) 2 (tfree) 2 



a = 2M 



3(ffree + tpause) 

and intercept 

„ 2t free., (l*ee) 2 (tfree) 3 
p = M = -. 

3 9(tf re e + tpausej 



(7) 



(8) 



(9) 



This is an important difference to other models of 
mean square displacement such as Fiirth's equation, 



which is frequently applied to cell migration data 
[16,17]. This equation reads 



E[D 2 (t)] = 2M{t - P{1 - e- t/p )) 



(10) 



where P is a parameter called the persistence length. 
Hence, in the Beauchemin Model, the persistence of 
migration carries on only for the predefined time £ pause + 
tf ree , after which it is immediately and completely lost. On 
the other hand, persistence in Fiirth's equation gradually 
decays over time following an exponential term. This is 
why we previously used Fiirth's equation to estimate the 
cell motility coefficient from short-term migration data 
[7]. However, whether this really leads to more accurate 
results remains to be determined. 

A second notable consequence of the above result is 
that although the Beauchemin model has three micro- 
scopic parameters, its step-wise variance has only two 
degrees of freedom. Consequently, while the value of 
£ P ause would become apparent when analyzing the track 
of a single simulated cell at a very high resolution, it 
does no longer matter when a large population of cells 
that migrate asynchronously is analyzed, as shown by 
Proposition 12: Given a fixed tf lee , there are infinitely 
many combinations of £ pause and v free that yield an iden- 
tical motility coefficient and thereby an identical 
expected square displacement of the population. For 
example, the empirical mean square displacement of the 
best-fitting triplet (0.5, 2.0, 18.8) is identical (up to some 



r \ 

exact variance 

at - /3 




time 

Figure 5 Forming the population average. Expected mean square displacement of observed particle ensembles as of Proposition 12. The 
pulses of the individual variances (see Figure 4) average out, and the exact expression converges quickly to a linear form with slope a = 2M 
and intercept /j = 2Mf free /3. One degree of freedom from the individual variance is lost due to the averaging process. 
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rounding) to that of the triplet (9.32, 2.0, 40). Hence, we 
cannot expect to get reliable estimates for all three para- 
meters tfrea. £ P ause and Vf ree from fitting the model to 
mean square displacement data. This is consistent with 
the observation by Beauchemin et al. that the parameter 
fpause can always be set to 0 without much influence on 
the quality of the best fit. 

However, Beauchemin et al. fitted their simulations to 
the mean displacement data, and the above result does 
not imply that the mean displacement too has just two 
degrees of freedom. Perhaps surprisingly, simulations 
show that this is in fact not the case (Figure 6): While 
mean square displacement plots generated with using the 
above parameter triplets and our analytical solutions are 
indeed consistent with each other, the mean displace- 
ment plots are significantly different. Hence, the mean 
displacement plot appears to preserve more information 
on the underlying migration process than does the mean 
square displacement plot, at least when applied to cell 
tracks simulated using the Beauchemin model. 

Confinement ratio 

The confinement ratio is a measure of track straight- 
ness. In the context of lymphocyte migration analysis, 
it is also called "chemotactic index" or "meandering 
index" [2,3]. The confinement ratio is usually defined 
as the quotient of a path's displacement from the ori- 
gin versus the complete length of the path. If the cell 
changes directions very often, the confinement ratio 
will tend to zero over time, while for directed 



migration without a random component, the confine- 
ment ratio is equal to 1. 

Instead of the confinement ratio, we consider here the 
squared confinement ratio, which can be easily determined 
for the Beauchemin model from the previous result. 

Definition 13 (Squared confinement ratio). Consider a 
particle moving according to the Beauchemin model and 
let D (t) denote the squared distance between the parti- 
cle position at time t and its starting point. Then we 
define the squared confinement ratio C (t) as follows: 



c 2 {t) 



DHt) 



By assuming i pause = 0 and applying the results from 
the previous section, we can determine the expected 
confinement ratio of the Beauchemin model as follows: 



E [C 2 (t)l - 



t 2 xL 



2Mt - 2Mt fce x 



(fee _ tfee 
it it 2 

it 9 X 

tta _ 1 
it 9 * 



§(&) 3 -(i) 2 +(£) — 



t > (free 
; - 3 + % t < tftee 

t > (fee 



1 



3 max(t, tfr ee ) 



- 3 t < W 



3 max (t, I fee ) 2 




Vthnc (V min) -\Aimc (Vmin) 

Figure 6 Mean displacement can be more informative than mean square displacement. The pause parameter of the Beauchemin model is 

not essential for the mean square displacement, but for the mean displacement. Left: Root mean square displacement curves generated for two 

populations with f free = 2 min and t pause = 0.5 min, v free = 18.8 /jm/min (red squares and error bars) or f pause = 9.32 min, i/ free = 40 ^m/min (blue 

circles and error bars). Simulations (N = 10000) with these two parameter yield the same results, as predicted by Proposition 12 (black line). Error 

bars show the standard error of the mean. Right: Although there is no difference with respect to the mean square displacement, the mean 

displacement curves of the two simulated populations differ. 
* J 
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Note that E [C 2 (0)] = |, because the path length tVf ree 
is measured in three dimensions while the expected 
observed square distance is projected to one dimension. 

Beltman et al. [3] have argued that the confinement 
ratio is not well suited for comparisons between differ- 
ent experiments because it tends to zero over time, and 
is therefore sensitive to the duration of the experiment 
from which it was estimated. They suggested to instead 
normalize the confinement ratio by multiplying it with 
the square root of time, such that it converges to a con- 
stant value. Applying this normalization to the squared 
confinement ratio of the Beauchemin model, we obtain 

(xE[C 2 (t)] = \ -ft ftee 1 * ^ . V (11) 

1 V " 3 max(t, t free ) V 3 max (t, £ free ) 2 / K ' 

which for t > t free simplifies to 
t x E [C 2 (£)] = tfree(3t ~ 1) . (12) 

Thus, the normalized squared confinement ratio con- 
verges to if r ee/3. When there is evidence that the 
observed migration process really behaves roughly like 
the Beauchemin model, this could be a good rule of 
thumb for quick estimations of the persistence time 
from experimental data. 

Simulating biased migration 

For our recent study of the detection limits of two- 
photon imaging [7], we extended the Beauchemin model 
for simulating three different types of biased migration, 
for which we used the term taxis modes [18] to indicate 
that such biased migration in cells usually occurs in 
response to an external stimulus. These three taxis 
modes were orthotaxis, where the movement speed is 
adjusted, klinotaxis, where the duration of the "free runs" 
is changed, and topotaxis, where the turning angle distri- 
bution of the simulated cell is changed. Thus, these mod- 
ifications can be viewed as semi-mechanistic ones that 
attempt to describe ways in which a cell could re-adjust 
its migration behaviour from pure random migration to 
biased random migration in response to an external sti- 
mulus. Before we review these modifications, we start by 
presenting a fourth method to simulate biased migration 
which is not mechanistically motivated, but is much 
easier to analyze. 

In all cases, taxis is spatially and temporarily uniform 
and controlled by two parameters: The vector ^ g r3, 
assumed to have unitary length, sets the direction of the 
taxis, and the parameter p & [0, 1] determines the mag- 
nitude of taxis, with p = 0 equivalent to an unbiased 
random walk and p = 1 to a maximally biased random 
walk; the precise meaning of the parameter p will be 
defined for each case below. 



Simple phenomenological model 

A simple ad-hoc method to extend the Beauchemin model 
for simulating biased migration is the following: During 
the pause phase, the cell is no longer kept stationary, but 
moves into the bias direction ^ with speed p ■ v {lee . This 
model is very simple to analyze, because the modification 
is deterministic and so changes only the mean, but not the 
variance of each step. Hence, the motility coefficient 
remains the same. 

Proposition 14 (Taxis speed for the simple model). 
The length of the convection coefficient, \\C\\, of the sim- 
ple biased migration model is given by 

_ P ' v free • tpquse 
tfree + tpause 

Proof. During each pause phase, the cell moves a dis- 
tance p ■ v free • tpause into the bias direction. 1 1 C| | is then 
simply the product of this distance and the total fraction 
of time that the cell spends in the pause phase. □ 

There might be some use for the simple biased migra- 
tion model in situations where it is necessary to study the 
effect of bias independently of the motility coefficient. 
However, it is difficult to imagine how such a situation 
could arise biologically. Below we therefore discuss model 
modifications that have a more realistic motivation. For 
each of these modifications, it is possible to analytically 
derive the convection coefficient, i.e., the expected speed 
of the bias. 

Throughout the rest of this section, we let the vector ^ 
denote the orientation vector of the Beauchemin model 
(i.e., a vector randomly sampled from the unit sphere), 
and i, the bias direction. Moreover, we let (f i/, f ') 
denote the vector of random variables representing a step 
of the modified model, while (£ v, 0 denotes a step in the 
original model. 

Orthotaxis model 

In the orthotaxis model, the migration speed is no 
longer constant but depends on the chosen target direc- 
tion. The more consistent this direction is with the tar- 
get direction, the faster the simulated cell will move. 
Such bias could occur for instance in response to some 
external dragging force, and might be illustrated by ima- 
gining a random stroll through a city on a very windy 
day. 

Formally, we define the adjusted speed as follows: 

"f ree = Vfree (l +P (&, 3)) 

Here (•>•) denotes the scalar product and p the bias 
strength parameter. Note that the movement perpendi- 
cular to fjis not affected at all, while movement against 
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the direction jj is slower and in direction it is faster. 
More precisely, the relation between p and the bias 
speed v taxis is given by the following proposition. 

Proposition 15 (Speed of orthotaxis [7]). The convec- 
tion coefficient \ \C\ \ of the orthotaxis model is given by 



\\Q\ = v fm 



P tfr a 



Proof. We can assume without loss of generality that 
the bias is along the first dimension, i.e. 2 = (1, 0, 0)> 
and thus (jj, d) = f • We are interested in the expecta- 
tion E [£']; it is easy to see that the expectations E [t/] 
and E [f '] are still 0. Then we get 

f = ffreel = ffree(l + PH)H 

= l*ee(§ +PH 2 )- 
Because f is uniformly distributed on [-1, 1] (see 
proof of Proposition 10), the expected location E [§'] at 
the end of the step is 

E[f] = E[y free - tfee^+Pf 2 )] 

= V free • tf ree (E [|] + pE [? 2 ]) 

= "free • tfree • p/3. 

Dividing this expected location E [£'] by the step 
duration f free + i pauS e yields the claimed expression for 
HC||. □ 

Topotaxis model 

In the case of topotaxis, the cell adjusts its turning 
behaviour but keeps the speed constant. This also affects 
both the variance and the mean of the cell step distribu- 
tion. Topotaxis could be a likely bias mechanism for 
cells migrating on some sort of structural cell network 
in the tissue, like it was suggested to be the case for 
T cells in secondary lymphoid organs [19]. 

We implement topotaxis by skewing the distribution 
from which we pick the cell's orientation d = (f , v, £)• 
Recall that this is a uniform distribution on a sphere in 
the unbiased case. Such a distribution can also be gener- 
ated by sampling f uniformly from [-1, 1] and then 
picking v, f uniformly from a circle with radius 
y/l — £ 2 . Thus, the probability density function f^(x) is 
defined as follows: 

f t (x) = \ , X€ [-1,1]. 

For the modified model version, we define a parame- 
terized version of fa (x) as follows: 

1 + px 

/{'(*) = — T~. *e[-l,l]. 



Hence, for p = 0 we obtain the distribution 
f^i(x) = T , of the unbiased case, while for p = 0.5 we 
get the slightly skewed distribution (.r)= --^ and for 
p = 1 we obtain the maximally skewed distribution 

This modification of the Beauchemin model yields the 
same convection coefficient as in the previous case. 

Proposition 16 (Speed of topotaxis [7]). The speed \ \ 
C\ | of the directional motion induced by the topotaxis 
model is given by the equation 



\\C\\ = 



P tfre 



3(tjree + tpause) 



Proof. This is simply shown by calculating E [£'] for 
the probability density function defined above. a 

Klinotaxis model 

Lastly, we can also induce bias by modifying the dura- 
tion f free of the free run depending on the picked orien- 
tation 2- A well-known biological example of such bias 
is the "run-and-tumble" motion by which E. Coli 
searches for food: Because the bacterium is too small to 
sense concentration gradients in situ, it is capable of 
comparing concentration measurements along a persis- 
tent path until it figures out whether it is going in the 
right direction. When it senses that it is moving in the 
right direction (run), it will try to keep going, while 
otherwise it will stop its motion soon (tumble) and ori- 
ent itself in a new random direction. 

We implement klinotaxis by modifying the free run 
duration t {ree depending on the angle between 2 and ^ 
as follows: 

4ee = £ free(l +P(b, 2» 

Again, this modification leads to exactly the same bias 
speed as found for the previous two models. 

Proposition 17 (Speed of klinotaxis [7]). The speed 
1 1 C\ | of the directional motion induced by the klinotactic 
modification of the Beauchemin model is given by the 
equation 



\\Q\ = Vfre, 



P tfre, 



Proof. The proof follows along the lines of a similar 
model of E. Coli run-and-tumble motion discussed by 
de Gennes [20]. We start again by analyzing a single 
free run of a cell projected onto the direction of bias, 
which is assumed without loss of generality to be along 
the first dimension £ Again we need to determine E [£']. 
We start by writing f as 

%' = Vfree ■ I • tfree • (1 +P£) • 
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Now let us assume for a moment that £ pause = 0. It is 
easy to verify that E [igj = tfree. Thus, the average 
speed towards the bias direction over many successive 
steps is E [£'/ £ free]- A simular calculation as in the proof 
of Proposition 15 yields 



E[£7k. 



p/3. 



Now we arrive at the claimed expression by multiply- 
ing the above with the overall fraction of time that the 
particle spends in the free run phase, which is 

ffree/(ffree + fpause)- 



Note that the klinotaxis model no longer fits within 
our mathematical framework because the duration of a 



step in the model is no longer constant. Thus, the cen- 
tral limit theorem we applied in that section no longer 
applies, and a different central theorem would be 
needed to show formally that the model converges to a 
Brownian motion. 

Moreover, we point out that for all modifications 
discussed in this section except for the simple phe- 
nomenological one, the random motility component 
changes. Therefore, these modifications cause as a 
side-effect a modified motility coefficient. In general, 
the random motility component are even affected in 
a non-isotropic fashion, such that the motility coeffi- 
cient would need to be described using a 3 x 3 
matrix instead of a single scalar value. In Figure 7, 
we show some simulation results that illustrate this 
point. 
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Figure 7 Taxis skews the random motility component. Displacement of in silico cells when applying the four different methods of simulating 
biased migration along the X axis using the methods discussed in this paper, i.e., the simple phenomenological model (black), orthotaxis (red), 
topotaxis (blue), and klinotaxis (green). For each panel, N = 1000 cells were simulated using the parameters f pause = 0.5 min, r free = 2.0 min, and 
\/f r ee = 18.8 (jm/min. Cell displacements are shown after one hour of simulation time, and the taxis parameters p are in each case set such that 
the expected displacement is 300 pm per hour. The length of the scale bars is proportional to the dimension-wise variances of the cell 
displacements, demonstrating that these variances change when biased migration is introduced. 
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Conclusions 

In this article, we have compiled both old and new analyti- 
cal results on the Beauchemin model of lymphocyte 
migration, which was originally a model of purely random 
motility, but is now able to accommodate partially biased 
migration as well. An important practical consequence of 
these results is that validation of the model against experi- 
mentally determined motility coefficient, mean square dis- 
placement data, and confinement ratio curves does no 
longer need to be performed by simulation, because these 
quantities can be calculated directly from the model para- 
meters. This also allows for the use of standard numerical 
fitting procedures, such as implemented in popular com- 
puter algebra packages, to estimate the model parameters 
from such data by least-squares-fitting. An important 
observation was that when simulating a large population 
of unsynchronized cells, the model is equivalent to the 
ideal chain model of polymer physics, which is however 
only true with respect to the quantities named above and 
not for others like the mean displacement. Interestingly, 
this shows that we can sometimes infer information from 
mean displacement data that we cannot infer from mean 
square displacement data. However, an analytical deriva- 
tion of the mean displacement of the model has not yet 
been achieved, and can be expected to be complicated jud- 
ging from the experiences with similar models [21]. For 
now, we leave this as an important open problem. 

Subtle biased migration of lymphocytes can point to 
important biological functions or phenomena [22,23], 
and the use of mathematical models has helped to under- 
stand the quantitative impact of such effects [23,7] . The 
modifications discussed above facilitate simulations of 
such subtle biases with pre-defined speed, again remov- 
ing the need for fitting the bias strength parameter to 
data by using simulations. However, all but one of the 
discussed modifications also affect the random migration 
component, which is likely also true in biological scenar- 
ios where randomly migrating cells respond to an exter- 
nal stimulus. Therefore, another important task that 
remains to be solved is to analytically derive the motility 
coefficient matrices for the topotaxis, klinotaxis, and 
orthotaxis models. 
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